#ifndef BENCHMARK_RUNNER_HPP
#define BENCHMARK_RUNNER_HPP

#include "solverSettings.hpp"

/* for quick tests */
struct benchmarkRunner {
	
	hms::Time time;
	
	hms::Traversal traversal;
	hms::EdgeTraversal edgeTraversal;
	hms::CellTraversal cellTraversal;
	Index blocksize;
	bool
		useFriction,
		useSourceFields,
		useMuscl,
		useExport,
		useAtomicWrite;
	scalar 
		maxCourant;
	fs::path exportPath;

	benchmarkRunner( hms::Time&& t, hms::SolverSettings&& s ) :
		time { std::move(t) },
		traversal { s.traversal },
		edgeTraversal { s.edgeTraversal },
		cellTraversal { s.cellTraversal },
		blocksize { s.blocksize },
		useFriction { s.useFriction },
		useSourceFields { s.useSourceFields },
		useMuscl { s.useMuscl },
		useExport { s.useExport },
		useAtomicWrite { s.useAtomicWrite },
		maxCourant { s.maxCourant },
		exportPath { std::move(s.exportPath) }
	{}

	template<typename MeshType>
	void operator()( MeshType& mesh ){
		using namespace hms;

		mesh.blocksize(blocksize);
		mesh.printInfo();

		SWESolver<MeshType> solver;

		// solver properties
		solver.useAtomicWrite(useAtomicWrite);
		solver.useFriction(useFriction);
		solver.useSourceFields(useSourceFields);
		solver.useMuscl(useMuscl);
		solver.limiterFunction(LimiterFunction::minMod);
		solver.maxCourant(maxCourant);

		// traversal
		solver.traversal(traversal);
		solver.edgeTraversal(edgeTraversal);
		solver.cellTraversal(cellTraversal);

		if(traversal == hms::Traversal::edges) {
			std::cout << "Edge traversal: "
				<< magic_enum::enum_name( solver.edgeTraversal() ) << std::endl;
		} else {
			std::cout << "Cell traversal: "
				<< magic_enum::enum_name( solver.cellTraversal()) << std::endl;
		}
		
		std::cout << "\nFriction is enabled: " << std::boolalpha << useFriction
			<< "\nSource fields are enabled: " << useSourceFields
			<< "\nUsing second order reconstruction: " << useMuscl
			<< "\n";
		if(useAtomicWrite) std::cout << "Atomic Write IS used!" << std::endl;
		else std::cout << "Atomic Write IS NOT used!" << std::endl;
		#ifdef _OPENMP
		std::cout << "Number of OpenMP threads is set to: " <<  omp_get_max_threads() << "\n";
		#endif

		// fields
		solver.v() = Field{ mesh, 2 };
		solver.d() = Field{ mesh, 1 };
		solver.z() = Field{ mesh, 1 };
		solver.v().values().colwise() = Array2s{ 0.0, 0.0 };
		solver.z().values() = 0.0;

		enum Case {
			 finn, dambreak2D
		};

		constexpr Case runCase { dambreak2D };

		if constexpr ( runCase == Case::finn ){
		
			solver.d().values() = 0.1;

			scalar zPlateau { 0.1 };

			mesh.forEachCell( [&](const auto& cell){
			
				scalar x {cell.centroid()[0]};
				scalar y {cell.centroid()[1]};
			
				if(x >= 0.0 && x < 10.0 && y >= 0.0 && y < 10.0) {
					solver.z().values()(0, cell.index()) = zPlateau;
				}
			});

			// Index
			// 	nx {mesh.nCellsX()},
			// 	ny {mesh.nCellsY()},
			// 	cellsPerRow {nx + 2*mesh.nGhostLayers()},
			// 	cellsPerCol {ny + 2*mesh.nGhostLayers()};

			solver.n() = ConstantField<Uniform, MeshType>{ mesh, ArrayNs<1>{0.025} };

			if (useSourceFields){
				solver.r() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
				solver.i() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
			}

			// boundary conditions
			const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
			const BoundaryPatch& top    { mesh.boundaryPatch("top") };
			const BoundaryPatch& right  { mesh.boundaryPatch("right") };
			const BoundaryPatch& left   { mesh.boundaryPatch("left") };

			scalar vBoundary {0.5};
			solver.v().setBoundaryCondition( left,   OrthogonalFlowVelocity_velocity{}, vBoundary );
			solver.v().setBoundaryCondition( right,  FreeOutflow_velocity{} );
			solver.v().setBoundaryCondition( top,    FreeOutflow_velocity{} );
			solver.v().setBoundaryCondition( bottom, OrthogonalFlowVelocity_velocity{}, vBoundary );

			solver.d().setBoundaryCondition( left,   OrthogonalFlowVelocity_waterdepth{}, solver.v(), vBoundary );
			solver.d().setBoundaryCondition( right,  FreeOutflow_waterdepth{} );
			solver.d().setBoundaryCondition( top,    FreeOutflow_waterdepth{} );
			solver.d().setBoundaryCondition( bottom, OrthogonalFlowVelocity_waterdepth{}, solver.v(), vBoundary);

			Array1Xs bottomElevationLeft   { solver.z().domainValues(left) };
			Array1Xs bottomElevationBottom { solver.z().domainValues(bottom) };

			solver.z().setBoundaryCondition( left,   StaticValue{},bottomElevationLeft, false );
			solver.z().setBoundaryCondition( right,  StaticValue{}, Array1Ns<1>{0}, true );
			solver.z().setBoundaryCondition( bottom, StaticValue{},bottomElevationBottom, false );
			solver.z().setBoundaryCondition( top,    StaticValue{}, Array1Ns<1>{}, true );
			// solver.z().setBoundaryCondition( top,          StaticValue{}, bottomElevationTop, false );
			// solver.z().setBoundaryCondition( bottomWall,   StaticValue{}, bottomElevationWall, false );
			// solver.z().setBoundaryCondition( bottomInflow, StaticValue{}, bottomElevationInflow, false );
		} else if constexpr ( runCase == Case::dambreak2D ){
		
			solver.d().values() = 0.0;

			scalar dPlateau { 0.1 };

			mesh.forEachCell( [&](const auto& cell){
			
				scalar x {cell.centroid()[0]};
				scalar y {cell.centroid()[1]};
			
				if(x >= 0.5 && x < 1.0 && y >= 0.5 && y < 1.0) {
					solver.d().values()(0, cell.index()) = dPlateau;
				}
			});

			// Index
			// 	nx {mesh.nCellsX()},
			// 	ny {mesh.nCellsY()},
			// 	cellsPerRow {nx + 2*mesh.nGhostLayers()},
			// 	cellsPerCol {ny + 2*mesh.nGhostLayers()};

			solver.n() = ConstantField<Uniform, MeshType>{ mesh, ArrayNs<1>{0.025} };

			if (useSourceFields){
				solver.r() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
				solver.i() = SourceField<Uniform, MeshType>{ mesh, ArrayNs<1>{0} };
			}

			// boundary conditions
			const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
			const BoundaryPatch& top    { mesh.boundaryPatch("top") };
			const BoundaryPatch& right  { mesh.boundaryPatch("right") };
			const BoundaryPatch& left   { mesh.boundaryPatch("left") };

			solver.v().setBoundaryCondition( left,   FreeOutflow_velocity{} );
			solver.v().setBoundaryCondition( right,  FreeOutflow_velocity{} );
			solver.v().setBoundaryCondition( top,    FreeOutflow_velocity{} );
			solver.v().setBoundaryCondition( bottom, FreeOutflow_velocity{} );

			solver.d().setBoundaryCondition( left,   FreeOutflow_waterdepth{} );
			solver.d().setBoundaryCondition( right,  FreeOutflow_waterdepth{} );
			solver.d().setBoundaryCondition( top,    FreeOutflow_waterdepth{} );
			solver.d().setBoundaryCondition( bottom, FreeOutflow_waterdepth{} );

			// Array1Xs bottomElevationLeft   { solver.z().domainValues(left) };
			// Array1Xs bottomElevationBottom { solver.z().domainValues(bottom) };

			// solver.z().setBoundaryCondition( left,   StaticValue{},bottomElevationLeft, false );
			// solver.z().setBoundaryCondition( right,  StaticValue{}, Array1Ns<1>{0}, true );
			// solver.z().setBoundaryCondition( bottom, StaticValue{},bottomElevationBottom, false );
			// solver.z().setBoundaryCondition( top,    StaticValue{}, Array1Ns<1>{}, true );
			solver.z().setBoundaryCondition( left,   ZeroGradient{} );
			solver.z().setBoundaryCondition( right,  ZeroGradient{} );
			solver.z().setBoundaryCondition( bottom, ZeroGradient{} );
			solver.z().setBoundaryCondition( top,    ZeroGradient{} );
		}

		// initialisation
		solver.initialise();

		// simulation time & export
		// Time time;

		// time.startAt(0.0);
		// time.endAt(simulationTime);
		// time.writeAtTimeIntervals(useWriteTimeInterval);
		// if ( useWriteTimeInterval ) time.writeTimeInterval(writeInterval);
		// else time.writeStepInterval(static_cast<hms::Index>(writeInterval));
		
		// measurement timer
		Timer t;
		// time marching loop
		std::cout << std::setprecision(15) << std::fixed;
	
		// create binding file
		std::unique_ptr<vtk::xml::XMLDocument> bindingFile {(useExport) ? vtk::createBindingDoc() : nullptr };

		while ( !time.atEnd() ){
			
			if( useExport && time.atWriteTime()) {
				vtk::toVtk(exportPath,mesh,solver,time.nWritten());
				vtk::updateBinding(bindingFile.get(), mesh, time.current(), time.nWritten(), exportPath);
			}
			solver.advance( time );
			if ( time.atWriteTime() )
				std::cout << "Time: "
					<< std::setprecision(3) << time.current() << "\n";
		}

		if( useExport ){
			vtk::toVtk(exportPath,mesh,solver,time.nWritten());
			vtk::updateBinding(bindingFile.get(), mesh, time.current(), time.nWritten(), exportPath);
		}
		double elapsed { t.elapsed() };


		// execution times
		double execTimeSum {0};
		for ( auto execTime : solver.execTimings() )
			execTimeSum += execTime;

		// print execution times for different parts of the computation
		std::cout
			<< "\nComputation time: " << elapsed << "s"
			<< "\nSum of execution times logged by solver: " << execTimeSum << "s"
			<< "\n\t" << 100 * solver.execTimings()[0] / execTimeSum
				<< "% for time step calculation"
				<< " (" << solver.execTimings()[0] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[1] / execTimeSum
				<< "% for fluxes (including slope)"
				<< " (" << solver.execTimings()[1] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[2] / execTimeSum
				<< "% for sources " << ( useFriction ?
					std::string{ "(including friction)" } :
					std::string{ "(friction disabled)" } )
				<< " (" << solver.execTimings()[2] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[3] / execTimeSum
				<< "% for new state"
				<< " (" << solver.execTimings()[3] << "s)"
			<< "\n\t" << 100 * solver.execTimings()[4] / execTimeSum
				<< "% for updating boundaries"
				<< " (" << solver.execTimings()[4] << "s)"
			<< "\n";
		
		#ifdef SOLVER_TIME_FLUXES
		for ( double& d : solver.execTimesFlux() )
			d /= omp_get_num_procs();
		
		double fluxTimeSum {0};
		for ( auto fluxTime : solver.execTimesFlux() )
			fluxTimeSum += fluxTime;
		std::cout
			<< "\nFlux calculation total: " << fluxTimeSum << "s"
			<< "\n\t" << 100 * solver.execTimesFlux()[0] / fluxTimeSum
				<< "% for getting cell indices"
				<< " (" << solver.execTimesFlux()[0] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[1] / fluxTimeSum
				<< "% for getting field values"
				<< " (" << solver.execTimesFlux()[1] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[2] / fluxTimeSum
				<< "% for handling dry cells"
				<< " (" << solver.execTimesFlux()[2] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[3] / fluxTimeSum
				<< "% for hydrostatic reconstruction"
				<< " (" << solver.execTimesFlux()[3] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[4] / fluxTimeSum
				<< "% for calculating wave speeds"
				<< " (" << solver.execTimesFlux()[4] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[5] / fluxTimeSum
				<< "% for calculating fluxes"
				<< " (" << solver.execTimesFlux()[5] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[6] / fluxTimeSum
				<< "% for slope source term"
				<< " (" << solver.execTimesFlux()[6] << "s)"
			<< "\n\t" << 100 * solver.execTimesFlux()[7] / fluxTimeSum
				<< "% for registering flux"
				<< " (" << solver.execTimesFlux()[7] << "s)"
			<< "\n";
		#endif

	}
}; // end of benchmark runner

#endif